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We present spin-independent and spin-spin interquark potentials for the charmonium and 
charmed-strange mesons, which are calculated in 2-1-1 flavor lattice QCD simulations using the 
PACS-CS gauge configurations generated at the lightest pion mass {Mt^ ~ 156(7) MeV) with a 
lattice cutoff of a~^ ~ 2.2 GeV and a spatial volume of (3 fm)^. For the charm quark, we use a rel¬ 
ativistic heavy quark (RHQ) action with fine tuned RHQ parameters, which closely reproduce both 
the experimental spin-averaged mass and hyper-fine splitting of the l^ charmonium. The interquark 
potential and the quark kinetic mass, both of which are key ingredients within the potential de¬ 
scription of heavy-heavy and heavy-light mesons, are determined from the equal-time Bethe-Salpeter 
(BS) amplitude. The charmonium potentials are obtained from the BS wave function of IS char- 
monia {rjc and J/tf mesons), while the charmed-strange potential are calculated from the Ds and 
D* heavy-light mesons. We then use resulting potentials and quark masses as purely theoretical 
inputs so as to solve the nonrelativistic Schrodinger equation for calculating accessible energy levels 
of charmonium and charmed-strange mesons without unknown parameters. The resultant spectra 
below the DD and DK thresholds excellently agree with well-established experimental data. 


I. INTRODUCTION 

The heavy-quark (Q)-antiquark (Q) potential is an 
important quantity to understand properties of the 
heavy quarkonium states. Because the dynamics of 
heavy quarks with masses much larger than the QCD 
scale (Aqcd) is well described within the framework of 
nonrelativistic quantum mechanics [1]. Indeed the con¬ 
stituent quark potential models with a QCD-motivated 
QQ potential have successfully predicted the heavy 
quarkonium spectra and its decay rates below open 
charm thresholds [2-4]. 

In such nonrelativistic potential (NRp) models, the 
conventional heavy quarkonium states such as charmo¬ 
nium and bottomonium are well understood to be quark- 
antiquark pair bound by the Coulombic potential induced 
by a perturbative one-gluon exchange that dominates in 
short range, plus linearly rising potential that describes 
the phenomenology of confining quark interactions at 
large distances [2]. This potential is called the Cornell 
potential and its functional form is given by 

+crr +Vb (1) 

3 r 

where as is the strong coupling constant, a denotes the 
string tension and Vq is the constant term associated with 
a self-energy contribution of the color sources. In ad¬ 
dition to the spin-independent potential, the NRp mod¬ 
els include spin-dependent interactions, which resolve the 
degeneracy among spin-multiplets. The spin-dependent 
potentials appear as relativistic corrections in powers of 
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the relative velocity of quarks, and their functional forms 
are also determined by perturbative one-gluon exchange 
as the Fermi-Breit type potential [5] . A more direct con¬ 
nection to QCD is established by the modern approach of 
effective field theory called potential nonrelativistic QCD 
(pNRQCD) [6]. 

We would like to stress here that the functional forms 
of the QQ potentials except at long distances are ba¬ 
sically deduced by the perturbative approach. Fur¬ 
thermore all of parameters needed in the NRp mod¬ 
els, including a constituent quark mass mg, are phe¬ 
nomenologically fixed to reproduce the experimental 
heavy quarkonium masses [3, 4]. The phenomenologi¬ 
cal spin-dependent potentials based on the perturbative 
method would have validity only at short distances and 
also in the vicinity of the heavy quark mass limit. This 
fact could cause large uncertainties in predictions for the 
higher-lying states of the heavy quarkonium in the NRp 
models. 

In addition, many of the charmonium-like mesons 
have been announced by B-factories at KEK and SLAC, 
which are primarily devoted to the physics of CP vio¬ 
lation, also by Charm factories at BEPC and CESR, 
and Tevatron at Fermilab. These newly discovered 
state above the open charm threshold could not be 
simply explained as conventional charmonium states in 
the constituent quark description [7]. Indeed, the ex¬ 
istence of the charged Z states including two charged 
bottomonium-like states, ^(,(10610) and Zb(10650) [8] in¬ 
dicates that the charmonium-like XYZ mesons are good 
candidates for non-standard quarkonium mesons such 
as hadronic molecular states, diquark-antidiquark bound 
states (tetraquark states), or hybrid mesons [9]. 

To discriminate between standard and nonstandard 
mesons in a zoo of the charmonium-like XYZ mesons, 
it is essential to investigate the validity of the potential 



2 


description of the heavy-heavy and heavy-light mesons 
directly from first principles of QCD. In this paper, we 
thus aim to provide the central QQ potentials (the spin- 
independent potential and the spin-spin potential), which 
are determined through the Bethe-Salpeter (BS) am¬ 
plitudes of pseudoscalar and vector mesons in dynami¬ 
cal lattice QCD simulations with almost physical quark 
masses. 

In lattice QCD, understanding the properties of QQ 
interactions is one of the great historic milestones. The 
Wilson loop has been originally introduced as a non-local 
order parameter in Z 2 gauge theory by Wegner [10]. Sub¬ 
sequently, Wilson generalized it with continuous gauge 
groups and related it to the static potential between in¬ 
finitely heavy-quark and antiquark in QCD so as to prove 
the quark confinement in the strong coupling limit [11]. 
The static QQ potential determined from Wilson loops 
have been precisely determined by lattice QCD in the 
past decades. The lattice QCD calculations within the 
Wilson loop formalism support a shape of the Cornell 
potential [12]. 

On the other hand, the spin-dependent QQ potentials 
regarded as the relativistic corrections to the static po¬ 
tential can be determined within the framework of pN- 
RQCD. Although earlier quenched studies [13, 14] and 
full QCD studies [15, 16] did not enable us to determine 
the functional forms of the spin-dependent terms due to 
large statistical errors, a full set of the spin-dependent 
terms {i.e. spin-spin, spin-orbit and tensor terms) have 
been successfully calculated in quenched QCD with high 
accuracy by using the multilevel algorithm [17, 18]. 

It is worth mentioning that the multilevel algorithm 
employed in Refs. [17, 18] is not easy to be implemented 
in dynamical lattice QCD simulations. Furthermore, the 
leading spin-spin potential determined at ©(l/mg) in 
quenched QCD gives an attractive interaction for the 
higher spin states in the hyperfine multiplet [17, 18]. This 
contradicts with the spin-spin term of the Fermi-Breit 
type potential, which is described by a repulsive contact 
interaction. Although one might think that the inverse 
of the charm quark mass would be far outside the valid¬ 
ity region of the 1 /toq expansion, this issue still remains 
even at the bottom quark mass. 

We develop the new method proposed in our previ¬ 
ous works [19-21] in order to obtain proper interquark 
potentials at finite quark masses, which are indispens¬ 
able for the potential description of the charmonium and 
charmed-strange mesons. The interquark potential and 
the quark kinetic mass, both of which are key ingredients 
within the potential description, can be defined by the 
equal-time and Coulomb gauge BS amplitude through 
an effective Schrodinger equation [19]. This new method 
enables us to determine the interquark potentials includ¬ 
ing spin-dependent terms at finite quark masses from first 
principles of QCD, and then fix all parameters needed in 
the NRp models. In our previous works with quenched 
lattice simulations [19, 21], we demonstrated that both 
spin-independent central potential and spin-spin poten¬ 


tial calculated in the BS amplitude method reproduce 
known results calculated within the Wilson loop formal¬ 
ism in the mq —^ 00 limit. We read off from our QQ 
potentials, which may encode all orders of the Ijmq ex¬ 
pansion, that the IfmQ expansion scheme may have the 
convergence behavior up to the bottom sector, while the 
charm sector is far outside the validity region for this 
expansion [21]. Furthermore, we found that the higher 
order corrections beyond the next-to-leading order are 
inevitably required for the repulsive feature of the total 
spin-spin potential even at the bottom sector [21]. In ad¬ 
dition, there is no restriction to extend the new method 
to dynamical calculation [20]. Hereafter we call the new 
method as BS amplitude method. 

Once one gets the reliable QQ potentials, which con¬ 
tain both the spin-dependent contributions as well as the 
spin-independent central one, we can easily verify how 
well the potential description is satisfied in the heavy- 
heavy and heavy-light meson systems through solving the 
nonrelativistic Schrodinger equation with purely theoret¬ 
ical inputs. If the potential description is valid, many 
physical observables such as mass spectrum of heavy- 
heavy and heavy-light mesons and their decay rates are 
easily accessible as is in the NRp models. 

In this paper, we extend our previous work [20] done 
in 2-1-1 flavor lattice QCD simulations using the PACS- 
CS gauge configurations [22] in order to investigate the 
validity of the potential description of the heavy-heavy 
and heavy-light mesons. The simulated pion mass (M^ « 
156(7) MeV) is close to the physical point, while the sim¬ 
ulated K meson mass as Mk ~ 554(2) MeV is about 10% 
heavier than the physical value. Although the strange 
quark is slightly off the physical point, the parameters 
of clover fermions for the strange quark are chosen to 
be equal to those of the strange sea quarks used in gauge 
field generation. For the charm quark, we employ the rel¬ 
ativistic heavy quark (RHQ) action that can control large 
discretization errors induced by large quark mass [23]. 
The RHQ parameters in the action were calibrated to re¬ 
produce the experimental spin-averaged mass and hyper¬ 
fine splitting of the IS' charmonium. 

We first concentrate on the heavy-heavy systems so 
as to calculate the charm quark mass and the charmo¬ 
nium potential from the BS amplitudes of IS charmonia 
{r]c and J/ip mesons). We reuse the data, which were 
previously published in Ref. [20], and then perform a 
more elaborate analysis proposed in Ref. [21]. New anal¬ 
ysis significantly reduces systematic uncertainties on the 
shape of the charmonium potential at short distances 
due to the usage of the highly improved Laplacian op¬ 
erator ^. Once the charmonium potential and the charm 
quark mass are precisely determined, we can numerically 


^ Note that the binding energy of the low-lying charmonium states, 
which we may consider to be nearly Coulombic bound states, are 
very sensitive to details of the short-range interaction. 
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solve the nonrelativistic Schrodinger equation with such 
theoretical inputs and without additional parameters. 

We then extend our research to the Dg heavy-light 
meson systems to extract the strange quark mass and 
the charmed-strange potential from the BS amplitudes 
of two lightest charmed-strange mesons {i.e. the Dg and 
Z?* heavy-light mesons). We will then discuss the valid¬ 
ity of the potential description on both charmonium and 
charmed-strange mesons. 

This paper is organized as follows. In Sec. II, we 
briefly describe the methodology to calculate the spin- 
independent and spin-dependent interquark potentials 
from the BS amplitude of heavy-heavy and heavy-light 
mesons in lattice QCD simulations. In Sec. Ill we give 
the details of parameters used in our Monte Carlo simu¬ 
lations, and then discuss the charmonium mass obtained 
from the standard lattice spectroscopy with two point 
correlation functions of mesons. In Sec. IV, we show nu¬ 
merical results of the BS wave function, the quark kinetic 
mass mq, the spin-independent central and spin-spin po¬ 
tentials, calculated from dynamical lattice QCD simula¬ 
tions. In Sec. V, we show the charmonium mass spectrum 
obtained by solving the nonrelativistic Schrodinger equa¬ 
tion with the theoretical inputs determined from dynam¬ 
ical lattice QCD simulations at almost physical point, 
and finally discuss possible systematic uncertainties on 
the resulting energy spectrum of the charmonium states. 
In Sec. VI, we present the results from an application to 
the Dg heavy-light meson systems. In Sec. VII, we sum¬ 
marize and discuss all results and future perspectives. 

II. FORMALISM 

In this section, we will briefly review the BS ampli¬ 
tude method to calculate the interquark potential with 
the finite quark mass. This is an application based on 
the approach originally used for the hadron-hadron po¬ 
tential, which is defined through the equal-time BS am¬ 
plitude [24-33]. More details of determination of the in¬ 
terquark potential are given in Ref. [21]. 

For simplicity, we here consider the case of the heavy 
quarkonium QQ. An extension to the heavy-light me¬ 
son made of two non-degenerate quarks is easy. In lat¬ 
tice simulations, we measure the following equal-time QQ 
BS amplitude in the Coulomb gauge for the quarkonium 
states [34, 35]: 

0r(r) = ^(OlQ(a;)rQ(a; + r)\QQ- J^^), (2) 

X 

where r is the relative coordinate between quark and 
antiquark at a certain time slice t. The operator T ap¬ 
peared in Eq. (2) represents the Dirac 7 metrics, which 
specifies the spin and the parity of meson operators. For 
instance, with 75 and 7 ^, one can form the pseudoscalar 
(PS) and the vector (V) operators with = 0“ and 
= 1~, respectively. A summation over spatial co¬ 
ordinates X projects out corresponding states with zero 


total momentum. The r-dependent amplitude, is 

called BS wave function. The BS wave function can be 
extracted from the four-point correlation function 

Gr(r,t,ts) = ^ {0\Q{x,f)TQ{x + r,t) 

x,x' ,y' 

X [Q{x',tg)VQ{y',tg))^ JO) (3) 

at large time separation between the source (ts) and sink 
(t) locations (Jt — ts\/a 3> 1) [21]. Here, the gauge held 
configurations are necessarily fixed to the Coulomb gauge 
at both time slices t and ts. In the limit of r —)■ 0, the 
four-point correlation functions are reduced to the two- 
point correlation functions of mesons with a wall source. 
In this paper, we focus only on the S'-wave BS wave func¬ 
tion (rjc and J/tp for the charmonium and Dg and D* for 
the charmed-strange meson), obtained by an appropriate 
projection to the Af representation in cubic group [36]. 

Below the inelastic threshold the BS wave function 
satisfies an effective Schrodinger equation with a nonlocal 
and energy-independent interquark potential U [24, 37, 
38] 

- ^'?i’r(T’)-f j dr'U{r,r')(j)T{r') = Ercprir), (4) 

where /i is the reduced mass of the QQ system. The en¬ 
ergy eigenvalue Et of the stationary Schrodinger equa¬ 
tion is supposed to be Mr —2 toq. If the relative quark ve¬ 
locity V = j V/mgj is small as u <C 1, the nonlocal poten¬ 
tial U can generally expand in terms of the velocity v as 
U{r',r) = {V(r)+Vs(r)5Q-5g+VT(r)^i2+VLs(r)i-5+ 
0{v'^)}S{r'-r) where 512 = (Sq ■ f){SQ-r) - Sq ■ Sq/3 
with r = r/r, S = Sq + Sq and L = r x (—iV) [24]. 
Here, V, Vs) and Vls represent the spin-independent 
central, spin-spin, tensor and spin-orbit potentials, re¬ 
spectively. 

The Schrodinger equation for 5-wave states is simpli¬ 
fied as 

+ Vir) + Sq ■ 5gVs(r)| Mr) = ErMr) (5) 

at the leading order of the u-expansion. Here, we es¬ 
sentially follow the NRp models, where the J/ip state 
is purely composed of the 15 wave function. However, 
within this method, this assumption can be verified by 
evaluating the size of a mixing between 15 and ID wave 
functions in principle. 

The spin operator Sq ■ Sq can be easily replaced by 
its expectation values: —3/4 and 1/4 for the PS and V 
channels, respectively. Then, the spin-independent and 


^ For the charmonium system, the inelastic threshold implies the 
DD threshold, while the DK threshold is a counterpart in the 
charmed-strange meson system. 
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spin-spin QQ potentials can be evaluated through the III. LATTICE SETUP AND HEAVY 

following linear combinations of Eq.(5): QUARKONIUM MASS 


V{r) 


Vs{r) 


(3 V^(/)v(r) 1 

\4 (/)v(r) 4 (^ps(r) 

1 f V^.^v(r) V^0Ps(r) ] 

mg I (/'ps(»') i ’ 


E -I_ 

-^ave “ 

mq 



where Save = Afave — 2mQ and Ehyp = My — Mps. The 
mass Mave denotes the spin-averaged mass as \Mps + 
\My. The derivative is defined by the discrete Lapla- 
cian on the lattice. 

In the BS amplitude method, there is a room for opti¬ 
mizing the differential operator since the discrete Lapla- 
cian is itself build in the definition of the interquark po¬ 
tential. In Ref. [21], we showed that the discrete Lapla- 
cian operator defined in the discrete polar coordinates 
called r-Laplacian is more suitable than the naive one de¬ 
fined in the Cartesian coordinates from the viewpoint of 
the reduction of the discretization artifacts on the short- 
range behavior of the interquark potential. The latter 
was adopted in our earlier works [19, 20], while we use 
the r-Laplacian throughout this paper. For details of the 
discrete Laplacian operators, we will explain in Sec. IV. 

The quark kinetic mass is also an important quantity 
in the determination of the interquark potentials since 
Eq. (6) and Eq. (7) requires information of the quark ki¬ 
netic mass mg. In Ref. [19], we propose to calculate the 
quark kinetic mass through the large-distance behavior 
of the difference of “quantum kinetic energies” (the sec¬ 
ond derivative of the BS wave function normalized by the 
BS wave function) between the spin-singlet and -triplet 
states in the hyperfine multiplet. The most simple choice 
is of course a pair of and states. Contributions of 
the long-range confining force are canceled out in the dif¬ 
ference of “quantum kinetic energies”. Under a simple, 
but reasonable assumption as limr_>.oo Vs(r) = 0 which 
implies there is no long-range correlation and no irrele¬ 
vant constant term in the spin-spin potential, one may 
expect that the difference of “quantum kinetic energies” 
at long distances stems only from the hyperfine splitting 
energy E^yp. Therefore, the quark kinetic mass can be 
read off in the following way: 


A. 2 + 1 flavor PACS-CS dynamical gauge ensemble 

The computation of the interquark potential for the 
charmonium (cc) and charmed-strange (cs) system is car¬ 
ried out on a lattice x Nt = 32^ x 64 using the 
2 + 1 flavor PACS-CS gauge configurations [22]. The 
gauge fields are generated by non-perturbatively 0{a)- 
improved Wilson quark action with csw = 1-715 [39] 
and Iwasaki gauge action at fi = 1.90 [40], which cor¬ 
responds to a lattice cutoff of a~^ = 2.176(31) GeV 
(a = 0.0907(13)fm) [22]. The spatial lattice size is of 
about Nga ~ 3 fm. The hopping parameters for the light 
sea quarks {K„d,Ks}={0.13781, 0.13640} give a pion mass 
of Mtt = 156(7) MeV and a kaon mass of Mk = 554(2) 
MeV [22]. Simulation parameters of dynamical QCD 
simulations used in this work is summarized in Table I. 
Although the light sea quark masses are slightly off the 
physical point, the systematic uncertainty due to this fact 
could be extremely small in this project. Our results are 
analyzed on all 198 gauge configurations, which are avail¬ 
able through International Lattice Data Grid and Japan 
Lattice Data Grid Gauge configurations is fixed to the 
Coulomb gauge. 


B. Relativistic charm quark 


In order to control discretization errors induced by 
large charm quark mass, we employ the relativistic heavy 
quark (RHQ) action [23] that removes main errors of 
0{\p\a), ©((tooo)") and 0{\p\a{mQa)'^) from the on-shell 
Green’s function. The RHQ action is the anisotropic ver¬ 
sion of the 0{a) improved Wilson action with five param¬ 
eters Kc, V, Ts, cb and ce, called RHQ parameters (for 
more details see Ref. [23, 41]): 

-S'rhq = ^ Qix) (moa + joDo + v‘y ■ D 

X 

-^a{DY 

+ ^ ^csaaijE^j + ^ l^CEaaQiFQ^Q{x) (9) 

ij * 


mq 


lim - 

r-->oo Ehyp 


(V^Mr) 

I <(>v(?') 


(l^psir) J 


( 8 ) 


The idea has been numerically tested, and the assump¬ 
tion of limr_>.oo Vs(r) = 0 is indeed appropriate in 
QCD [19]. We thus estimate the quark kinetic mass from 
asymptotic behavior of the right-hand side of Eq. (8) in 
long-distance region. 


where the Wilson parameter for the time derivative is 
set to be rj = 1 and the bare quark mass is related to 

the hopping parameter as amg = ^ - rt — Srg. The 

RHQ action utilized here is a variant of the Fermilab 
approach [42] (See also Ref. [43]). 


® International Lattice Data Grid/Japan Lattice Data Grid, 
http://www.jldg.org. 
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TABLE I: Parameters of 2 + 1-flavor dynamical QCD gauge field configurations generated by the PACS-CS collaboration [22]. 
The columns list number of flavors, lattice volume, the /3 value, hopping parameters for light and strange quarks, approximate 
lattice spacing (lattice cut-off), spatial physical volume, pion mass, and number of configurations to be analyzed. 


Nf 

Af X Nt 

P 


Ks 

a [fm] (a ^ [GeV]) 

NaU [fm] 

M.^ [MeV] 

# configs. 

2-bl 

32^ X 64 

1.9 

0.13781 

0.13640 

0.0907(13) (« 2.18) 

2.90(4) 

«156 

198 


TABLE II: The hopping parameter kq and RHQ parameters 
{v, Va, cb and ce) used for the charm quark. 



V 

Vs 

Cb 

Ce 

0.10819 

1.2153 

1.2131 

2.0268 

1.7911 


The parameters rg, cb and c^; in RHQ action are de¬ 
termined by tadpole improved one-loop perturbation the¬ 
ory [41]. For u, we use a nonperturbatively determined 
value, which is tuned by reproducing the effective speed 
of light Ceff to be unity in the dispersion relation E'^ (p^) = 
-I- CgffP^ for the spin-averaged 15'-charmonium state, 
since the parameter v is sensitive to the size of hyper- 
fine splitting energy [44]. We choose the value of Kc to 
reproduce the experimental spin-averaged mass of l^- 
charmonium states M^^p{1S) = 3.0678(3) GeV. To cal¬ 
ibrate RHQ parameters, we employ a gauge invariant 
Gauss smearing source for the standard two-point cor¬ 
relation function with four finite momenta. As a result, 
the relevant speed of light in a energy-momentum disper¬ 
sion relation -|- c^gP^ is consistent with unity 

within statistical uncertainties: c^g = 1.04(5) for the 
spin-averaged state [20]. Our chosen RHQ parameters 
are summarized in Table H. 

Using tuned RHQ parameters, we compute the two va¬ 
lence quark propagators with wall sources located at dif¬ 
ferent time slices tg/a = 6 and 57 to increase statistics. 
Two sets of two- and four-point correlation functions are 
constructed from the corresponding T operator with the 
charm quark propagator, and folded together to create 
the single correlation function. Dirichlet boundary con¬ 
ditions are imposed for the time direction at t/a = 0 
and 63 to eliminate unwanted contributions across time 
boundaries. 


C. Charmonium spectroscopy from two-point 
functions 
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Fig. 1 shows the effective mass of the S'-wave {rj and 
J/i/j) and P-wave {xco, Xci and he) charmonium states 
calculated from the dynamical lattice QGD simulation. 
These five charmonium states are classihed with quantum 
numbers and corresponding operators F as shown in 
Table HI. A effective mass is defined as 


FIG. 1: Effective mass plots for rje (upper panel), J/i/i (center 
pannel) and IP charmonium states (XcO, Xci and he) (lower 
pannel). Charmonium states are specified in legend. Solid 
lines indicate fit results and shaded bands display the fitting 
ranges and one standard deviations estimated by the jackknife 
method. 


Gr{t, ts) 

Gr{t + 1, ts) 


Mr{t) = log 


( 10 ) 
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TABLE III: Masses of low-lying charmonium states calcu¬ 
lated from two-point functions, the spin-averaged mass and 
hyperfine splitting energy of 15 charmonium states. Five 
charmonium states are classified with quantum numbers 
and corresponding operators F. The fitting ranges and values 
of x^/d.o.f. are also included. Results are given in units of 
GeV. 


state 

(J^^) 

r 

fit range 

mass [GeV] 

xVd.o.f. 

Vc 


75 

[33:47] 

2.9851(5) 

0.70 

J/^ 

(1-+) 

li 

[33:47] 

3.0985(11) 

0.62 

Mave(15) 

— 

— 

— 

3.0701(9) 

— 

^hyp(15) 

— 

— 

— 

0.1138(8) 

— 

XcO 


1 

[14:26] 

3.3928(59) 

0.66 

Xcl 

(1++) 

757i 

[14:26] 

3.4845(62) 

1.03 

he 

(1+-) 

7Gt 

[14:26] 

3.5059(62) 

0.63 


where Grit, tg) is the two-point function obtained by set¬ 
ting r to be zero in the four-point function Gr{r, t, tg) de¬ 
fined in Eq. (3). Each effective mass plot shows a reason¬ 
able plateau in the range 33 < t/a < 47 for S'-wave char¬ 
monium states and 14 < t/a < 26 for P-wave charmo¬ 
nium states. We estimate masses of the five charmonium 
states by a constant fit to the plateau over time ranges 
shown in Table III. A correlation among effective masses 
measured at various time slices is taken into account by 
using a covariance matrix in the fit. An inversion of the 
covariance matrix is performed once for average and it is 
used for each jackknife block. The statistical uncertain¬ 
ties indicated by shaded bands in Fig. 1 are estimated by 
the jackknife method. In Table II, we summarize resul¬ 
tant charmonium masses together with fit ranges used in 
the fits and values of per degrees of freedom (d.o.f.). 
Note that all masses calculated in this study are ob¬ 
tained from the Coulomb-gauge wall source propagator, 
while gauge-invariant Gaussian smeared source was used 
for results of charmonium masses compiled in Table I of 
Ref. [19]. 

Low-lying charmonium masses calculated below DD 
threshold are all close to the experimental values, though 
the hyperfine mass splitting Mhyp = 0.1124(9) GeV is 
slightly smaller than the experimental value, = 

0.1166(12) GeV [45]. The similar value of the hyper¬ 
fine mass splitting is reported even on the exact physical 
point in Ref. [44, 46]. Note that here we simply neglect 
the disconnected diagrams in all two-point correlation 
functions. The several numerical studies reported that 
the contributions of charm annihilation to the hyperfine 
splitting of the 15'-charmonium state is sufficiently small, 
which is of order 1-4 MeV. [47-49]. At the charm sec¬ 
tor, the effect of the disconnected diagrams on the char¬ 
monium, especially on the vector state, is perturbatively 
expected to be small due to Okubo-Zweig-Iizuka suppres¬ 
sion. 



FIG. 2: The reduced QQ BS wave functions of the r]c (circles) 
and J/ip (squares) states, shown as a function of the spatial 
distance r. The data points are taken along r vectors which 
are multiples of three directions (1,0, 0), (1,1, 0) and (1,1,1). 


IV. DETERMINATION OF INTERQUARK 
POTENTIAL 

A. QQ BS wave function 

We calculate the BS wave functions only for 5-wave 
states {r]c and J/ip). This is simply because the Goulomb- 
gauge wall source ^ adopted in this study is not suitable 
for studying the wave function of P-wave states, whose 
spatial part is odd under spatial reflection. 

Fig. 2 shows the QQ BS wave functions of 15- 
charmonium states {r]c and J/ip states). The BS wave 
functions are defined by Eq.(2) with a normalization con¬ 
dition oi'Yl^’T — 1- We use the reduced wave function 
Mr('r) for displaying the wave function: Mr('r) = r(j)Y{r). 
We practically take a time-average of the BS wave func¬ 
tion at fixed r over the range 33 < f/a < 47, where effec¬ 
tive mass plots for 15-charmonium states show excellent 
plateaus and excited state contaminations should be neg¬ 
ligible. To resolve the strong correlations between data 
of the BS wave function at different time slices, we take 
into account the covariance matrix during the averaging 
process over the time slice. 

We find that a breaking of rotational symmetry for 
the QQ BS wave functions is sufficiently small in our cal¬ 
culation. The resulting wave functions become isotropic 
with the help of a projection to the sector of the cubic 
group that corresponds to the 5-wave in the continuum 
theory (Fig. 2). All data points of the QQ BS wave func¬ 
tions calculated in the three different directions fall onto 
a single curve. 


^ Clearly, the spatial part of the meson operator constructed from 
a local quark bilinear operator with the wall source, where the 
quark operator is summed over all spatial sites at given time 
slice, belongs to the trivial A'^ irreducible representation of the 
cubic group. The plus sign in superscript indicates even parity. 
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FIG. 3: The determination of quark kinetic mass within 
the BS amplitude method. The values of — (V^c^v/i^v — 
V^0ps/0ps)/7?hyp as a function of the spatial distance r are 
shown in this figure. The quark kinetic mass rriQ is obtained 
from the long-distance asymptotic values of — (V^0v/'^v ~ 
^'^4‘Ps/4‘Ps)/Ehyp- Horizontal solid line indicates a value of 
quark kinetic mass obtained by fitting a asymptotic constant 
in the range 0.54 fm < r < 1.10 fm. A shaded band indicates 
a statistical error estimated by the jackknife method. 


The spatial lattice extent Nga « 2.9 fm is sufficiently 
large enough to study the 15'-charmonium system. In¬ 
deed, the BS wave functions shown in Fig. 2 are local¬ 
ized around the origin and vanished at r > 1.1 fm. It 
suggests that the QQ BS wave functions for the rjc and 
J/V' states fair enough fit into the box N^. Needless to 
say, the localized wave functions is interpreted as a sign 
of bound states. This fact however reminds us that the 
interquark potential can be deduced within the interior 
of the hadron due to its localized wave function. This is 
simply because that the signal-to-noise ratio in the cal¬ 
culation of of Eq. (6)-(8) is getting worse outside 

the spatial size of the hadron. 

Other important information can be read off from 
Fig. 2. The spatial size of the J/')/’ state is slightly larger 
than that of the rjc state. This indicates that there is 
a repulsive spin-spin interaction near the origin for the 
higher spin states. It is consistent with the pattern of 
level ordering for the hyperfine multiplet. The spin-spin 
charmonium potential will be discussed in more detail 
later. 


B. quark kinetic mass 

In our formalism, the kinetic mass of the charm quark 
is determined self-consistently within the BS amplitude 
method as well [19]. According to Eq. (8), the quark ki¬ 
netic mass can be evaluated from an asymptotic behav¬ 
ior of the quantity — (V^(/)v/</f'v ~ ^‘^(t>Ps/4>Ps)/Ehyp at 
long distances. For the discrete Laplacian operator V^, 
we use r-Laplacian, which is defined in polar coordinates 


as follows: 

r 2a 

^ (j)r{r -b d) -b (j)rir - d) - 2(j)r 

a? 

where r is the absolute value of the relative distance as 
r = jrj, and d is a spacing between grid points along 
differentiate directions. In the on-axis direction r oc 
(1,0,0) (labeled by “on-axis”), two off-axis directions 
r oc (1,1,0) (labeled by “off-axis I”) and r oc (1,1,1) 
(labeled by “off-axis 11”), the effective grid spacings cor¬ 
respond to d = a, y/2a and respectively. The dif¬ 
ference of ratios V^0r/?i>r at each r are obtained by a 
constant hts to the lattice data with reasonable y^/d.o.f. 
values over the range of time slices where two-point func¬ 
tions exhibit the plateau behavior (33 < t/a < 47). 

Fig. 3 illustrates the determination of quark kinetic 
mass tuq for the charmonium system. The value of 
mg can be determined from an asymptotic value of 
- V^(l)ps/4>Ps)/Ehyp in the range of 6 < 
r/a < 7V^ (0.54 fm < r < 1.10 fm), where V 5 (r) should 
vanish. In this study, three data sets are obtained from 
three directions: on-axis, off-axis I and off-axis II, are sep¬ 
arately analyzed so as to expose the size of the possible 
lattice discretization artifacts. On each data set, a value 
of niQ is obtained by a constant fit to a long-distance 
asymptotic value over the range as described above. Fi¬ 
nally we average them over three directions, and then 
obtain mg = 1.784(23)(6)(20) GeV. The first error is 
statistical, given by the jackknife analysis. In the second 
error, we quote a systematic uncertainty due to rota¬ 
tional symmetry breaking by taking the largest difference 
between average value and individual ones obtained for 
specific directions. The third ones are systematic uncer¬ 
tainties stemming from choice of tmin in tbe averaging 
process over the time-slice range tmin/a < t/a < 47. The 
minimum time-slice tmin/® is varied over range from 33 
to 41 and then take a largest difference from the preferred 
determination of mg. 

The charm quark mass obtained in this study is some¬ 
what heavier than the usual quark kinetic mass in the 
NRp models. For example, the quark kinetic mass 
adopted in Ref. [4] is about 17% smaller. This differ¬ 
ence however should not be taken seriously, because the 
value of mg in the NRp models highly depends on a con¬ 
stant term Vq of the Cornell potential, and Vq is actually 
forced to be zero in many of the NRp models. In ad¬ 
dition, the spatial profile of the spin-spin potential from 
lattice QCD is slightly different from the one used in the 
NRp models as we will discuss later. 

C. Spin-independent interquark potential 

Once the quark kinetic mass is determined, we can cal¬ 
culate the central spin-independent and spin-spin char¬ 
monium potentials from the QQ BS wave function 













0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 


r [fm] 

FIG. 4: Central spin-independent and spin-spin charmonium 
potentials calculated from the BS wave functions in the dy¬ 
namical QCD simulation with almost physical quark masses. 
In the upper panel, we show the spin-independent potential 
V{r). A solid (dot-dashed) curve is the fit results with the 
Cornell (Cornell plus log) form. The shaded bands show sta¬ 
tistical uncertainties in the fitting procedure where the jack¬ 
knife method is used. Note that the spin-averaged mass of 
IS'-charmonium states Save is not subtracted in this figure. 
A horizontal line indicates the level of open-charm (D^D^) 
threshold « 3729 MeV. In the lower panel, we show the spin- 
spin potential Vs(r). A solid (dot-dashed) curve corresponds 
to fitting results with exponential (Yukawa) form. The inset 
shows a magnified view. In both panels, the phenomenolog¬ 
ical potentials adopted in a NRp model [4] are also included 
as dashed curves for comparison. 


TABLE IV: Summary of the Cornell parameters and the 
quark mass determined by the BS amplitude method. 
For comparison, ones adopted in a phenomenological NRp 
model [4] and ones of the static potential obtained from 
Polyakov line correlations are also included. In the first col¬ 
umn, the quoted errors indicate the sum of the statistical and 
systematic added in quadrature. 



This work 

Polyakov lines 

NRp model [4] 

A 

0.713(83) 

0.476(81) 

0.7281 

[GeV] 

0.402(15) 

0.448(16) 

0.3775 

mq [GeV] 

1.784(31) 

OO 

1.4794 


through Eq. (6) and Eq. (7). Eirst, we show a result 
of the spin-independent charmonium potential V (r) in 
the upper panel of Fig. 4. The constant energy shift 
Save is not subtracted in this figure. As we reported 
in Ref [20], the charmonium potential calculated by the 
BS amplitude method from dynamical lattice QCD sim¬ 
ulations properly exhibits the linearly rising potential at 
large distances and the Coulomb-like potential at short 
distances. At first glance, the data points of the charmo¬ 
nium potential obtained from lattice QCD roughly follow 
the phenomenological potential used in the NRp mod¬ 
els, which is represented by the dashed curve. Neverthe¬ 
less, the data points at short distances are slightly off the 
dashed curve. In addition, a string breaking-like behavior 
is found in the range r <1.1 fm, where the charmonium 
potential reaches the level of open-charm threshold. We 
will discuss this point later, and for a while we concen¬ 
trate only on data points in the range r < 1.1 fm, where 
the linearly rising potential is clearly visible. 

For more close comparison, as a first step, we simply 
adopt the Cornell parametrization to fit the data of the 
spin-independent central potential: 

V{r) = - h err-I-Vb (12) 

r 

with the Coulombic coefficient A, the string tension u, 
and a constant Vq. All fits are performed individu¬ 
ally for each three directions over the range [rmi„/a : 
J'max/a] = [4 : 7^3], where r^ax « 1.1 fm. We min¬ 
imize the x^/d.o.f. with the covariance matrix and get 
the Cornell parameters of the charmonium potential as 
A = 0.713(26)(38)(31)(62) and = 0.402(6)(4)(9)(9) 
MeV with y^/d.o.f. « 3.2. The first error is statistical, 
and the second, third and forth ones are systematic un¬ 
certainties due to the choice of data points taken from 
three directions, and variations of t„iin and r^-m, respec¬ 
tively. 

The resulting Cornell parameters are summarized in 
Table IV. We also include both phenomenological ones 
adopted in a NRp model [4] and ones of the static poten¬ 
tial obtained from Polyakov-line correlators. The latter 
is estimated using the same method as in Ref. [22]. Ad¬ 
ditionally, we calculate the Sommer parameter defined 
as To = a/(1.65 — A)/a, and then obtain the value of 
ro = 0.476(6)(11)(3)(6) fm, which is fairly consistent 
with the value quoted in Ref. [22]. 

From our previous research in quenched QCD [21], the 
finite mg corrections could be encoded into the Cornell 
parameters. Indeed, as shown in Table IV, in the charmo¬ 
nium potential from the BS wave function, a Coulomb¬ 
like behavior is enhanced and the linearly rising force is 
slightly reduced due to finite charm quark mass effects 
in comparison to the conventional static potential from 
Wilson loops or Polyakov-line correlators. Furthermore, 
a gap for the Cornell parameters between the static and 
the phenomenological potentials seems to close by our 
new approach, which nonperturbatively accounts for a 
finite quark mass effect. 
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Here, we give a few technical remarks on the system¬ 
atic uncertainties on the Coulombic coefficient A, which 
highly depends on the choice of minimum value rj^in of 
the fitting window compared to the string tension a. This 
is simply because the linear part in the Cornell potential 
parametrization is dominated in the region in which we 
have data points. Indeed, as shown in the upper panel of 
Fig. 4, a solid curve, which corresponds to fitting results 
with the Cornell potential form, does not describe well 
data points outside the range [rmin : ’’max] used in the fit. 

In order to provide an adequate fit to the data points at 
shorter distances, we employed several alternative fitting 
forms. We found that a simple extension of the Cornell 
potential can describe the behavior of our charmonium 
potential reasonably well. An alternative htting form 
is given such that a log term is added to the Cornell 
potential: 

V(r) = - — + ar + Vo + B log(rA) (13) 

r 

where the value of A is simply set to be lattice cutoff a~^. 
Such a logarithmic r-dependence may appear in the lead¬ 
ing Xjmq correction to the static potential as reported in 
Ref. [50]. Moreover, as reported in Ref. [51], the charmo¬ 
nium potential obtained from the BS amplitude is consis¬ 
tent with the QQ potential obtained in the Wilson-loop 
approach within errors, when a leading Xjmq correction 
calculated in Ref. [52] is added to the static potential 
from Wilson loops. 

A fit with the “Cornell-plus-log” form (13) leads 
to the values of A = 0.194(137)(33)(36)(66), 

^ = 0.300(38)(19)(20)(21) GeV and B = 

0.390(113)(20)(39)(61) GeV with the slightly smaller 
value of x^/d.o.f. « 2.3. We here chose the fitting range 
to be [rmin/a : ’’max/a] = [3 : 7-\/3] and used a covariance 
matrix for taking into account the correlation among 
data points in the fit. The quoted errors have the same 
meaning as described above. 

We also plot the fit result with the Cornell-plus-log 
form, which is represented by a dot-dashed curve, in 
upper panel of Fig. 4. The shaded band displays one 
standard deviation estimated by the jackknife method. 
The short-distance behavior of the charmonium poten¬ 
tial is better described by the Cornell-plus-log form than 
the Cornell form (solid curve). If compared with the 
phenomenological potential of the NR models, the shape 
of the fitted curve of the Cornell-plus-log form at long- 
distances are much in agreement with the NR models 
though the string tension a becomes a slightly smaller 
value compared with the phenomenological one. In this 
context, the inclusion of the log term into the Cornell 
form gives only a minor modihcation at long-distances 
as far as the data is accessible in this study. 

Finally, we would like to comment on the string 
breaking-like behavior appeared in the range r > 1.1 fm. 
Although in principle, string breaking due to the pres¬ 
ence of dynamical quarks is likely to be observed, the 
observed feature in this study is suspicious and unreli- 


TABLE V: Results of fitted parameters for the spin-spin po¬ 
tential with the exponential and Yukawa forms. The quoted 
errors are statistical only. In the case of the spin-spin poten¬ 
tial, we use only the on-axis data. 


Functional form 

a 

P 

X^/d.o.f. 

Exponential 

2.15(7) GeV 

2.93(3) GeV 

2.0 

Yukawa 

0.815(27) 

1.97(3) GeV 

1.7 


able. As mentioned previously, the signal-to-noise ratio 
on the quantity of ^‘^(pr/'pr becomes worse rapidly as 
the spatial distance r increases because of the localized 
nature of the BS wave function (j)r{r). Moreover, the 
lattice data of the potential near the spatial boundaries 
are also sensitive to the possible distortion of its spatial 
profile as finite volume effects. Therefore, at least, cal¬ 
culations of the higher charmonium near the open charm 
threshold using a larger lattice is necessary for observing 
the string breaking in this sense. Their wave functions 
are extended until the string breaking sets in. 

We also emphasis that there might be another possi¬ 
ble reason for no evidence of string breaking from a view 
point of studies within the Wilson loop approach [53- 
56]. The string breaking in the static heavy quark po¬ 
tential can be observed only after inserting a opera¬ 
tor of light quark-antiquark to create the heavy-light 
meson-antimeson state (Qq){qQ), because the QQ cre¬ 
ation operator poorly overlaps with \{Qq){qQ)) state in 
Fock space [57-59] (See also Ref. [60] in the case of 
nonzero temperature). It is worth reminding that the 
static potential from Wilson loops is regarded as the “en¬ 
ergy eigenvalue” of the considering states. There would 
be nothing to change for the charmonium potential ex¬ 
tracted from the “stationary” wave function of the char¬ 
monium state, which is well defined in the BS amplitude 
method unless its energy level is above the open-charm 
threshold. 


D. Spin-Spin potential 

The lower panel in Fig. 4 shows the spin-spin charmo¬ 
nium potential obtained from the BS amplitude method 
with almost physical quark masses. The spin-spin poten¬ 
tial exhibits the short-range repulsive interaction, which 
is required to lead to energy gain for the higher spin 
state. Recall that the Wilson loop approach currently 
dose not achieve to reproduce the correct behavior of the 
spin-spin interaction. The leading-order spin-spin poten¬ 
tial classified in pNRQCD becomes attractive at short 
distances [17, 18]. Their calculation at next-to-leading 
order is unavailable at present. In contrast of the case 
of the spin-independent potential, the spin-spin poten¬ 
tial obtained from the BS wave function is absolutely 
different from a repulsive d-function potential generated 
by perturbative one-gluon exchange [5]. Such contact 
form oc 6{r) of the Fermi-Breit type potential is widely 
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adopted in the NRp models [3, 61]. 

The pointlike spin-spin interaction easily lifts the mass 
degeneracy between l^Pi state {he) and spin-weighted 
average of l^Pj states (xcj); M{l^Pj) = {M^^^ + 
- 1 - On the other hand, a finite-range 

interaction gives a non-zero, but small finite hyperfine 
splitting to the P- or higher-wave charmonia [7]. In the 
current experiments, however, the splitting Mhyp(lP) = 
M{l^Pj) — Mh^ for IP-charmonium states is not appre¬ 
ciably observed within experimental error. Here we quote 
-^hyp(l-P) = 0.02±0.19(stat) ±0.13(syst) MeV from the 
CLEO experiment [62, 63] (See also Ref. [64]). 

The QQ interaction is not entirely due to one-gluon 
exchange so that the spin-spin potential is not necessary 
to be a simple contact form oc S{r) [65-68]. This is shown 
to be true even for the ©(I/wq) spin-spin corrections in 
the Wilson-loop approach [17, 18], regardless of the sign 
issue. In the phenomenological side, the finite-range spin- 
spin potential described by the Gaussian form is adopted 
by some NRp model in Ref. [4], where many properties 
of conventional charmonium states at higher masses are 
predicted. This phenomenological spin-spin potential is 
also plotted in the lower plot of Fig. 4 for comparison. 
There is a slight difference at very short distances, al¬ 
though the range of spin-spin potential calculated from 
the BS amplitude method is similar to the phenomeno¬ 
logical one. 

To examine an appropriate functional form for the 
spin-spin potential, we try to fit the data with several 
functional forms, and explore which functional form can 
give a reasonable fit over the range of r/a from 2 to 7-\/3. 
As a results, the long-range screening observed in the 
spin-spin potential is accommodated by the exponential 
form and the Yukawa form: 

,, . \ _ / aexp(—/3r) : Exponential form , , 

( aexp(—/3r)/r : Yukawa form ^ ' 

All results of correlated fits are summarized in Ta¬ 
ble V. We also try to fit with the Gaussian form that 
is often employed in the NRp models, and it however 
gives an unreasonable y^/d.o.f. value. Note that we here 
use only the on-axis data which are expected to less suf¬ 
fer from the rotational symmetry breaking and the dis¬ 
cretization error, because fit results to the lattice data 
taken in each direction significantly disagree with each 
other [21]. We need the finer lattice to have a solid con¬ 
clusion to the shape of the spin-spin potential and the 
uncertainties due to the rotational symmetry breaking. 

V. CHARMONIUM MASS SPECTRUM FROM 
CHARMONIUM POTENTIAL 

Once the quark kinetic mass and the charmonium po¬ 
tentials are determined by first principles of QGD, we 
can solve the nonrelativistic Schrddinger equation defined 
with the theoretical inputs for the bound cc systems as 
same as calculations in the NRp models [2, 65, 66 , 69]. 


In the non-relativistic description, each charmonium 
state is generally labeled by a symbol ^'®+^[L]j, with the 
spin angular momentum (S' = 0,1, • • •), the orbital angu¬ 
lar momentum ([L] = S, P, D... corresponding to L = 
0,1, 2, • • •) and the total angular momentum {J — S(BL) 
quantum numbers. The notation is also used to 

classify the charmonium state. The parity (P) and the 
charge-conjugation {C) are given by P = (—1)'^+^ and 
C = (—within the non-relativistic description. 

Recall that all of the charmonium states below the 
open-charm threshold are experimentally well estab¬ 
lished [45]. The last missing IP-charmonium state, he, 
and also the first excited state of the pseudoscalar IS- 
charmonium state, r]e{2S), have already been observed 
in recent experiments [62, 63, 70-75]. 

In this section, we will discuss whether we can get the 
correct low-lying charmonium spectra within the hybrid 
approach between lattice QGD simulations and the NRp 
models in comparison to experimental data. In addition, 
we also perform a consistency check between two differ¬ 
ent methods in lattice QGD to verify the validity of our 
approach. One is of course the standard lattice spec¬ 
troscopy, where the mass information is extracted from 
the large-time asymptotic behavior of the two-point cor¬ 
relation functions, while another mainly uses the infor¬ 
mation about the spatial profile of the BS amplitudes. In 
this sense, these two methods are essentially different. 


A. Nonrelativistic Hamiltonian from lattice QCD 


For solving a nonrelativistic Schrddinger equation, the 
constant energy shift is irrelevant. We here introduce 
the energy-shifted potential in the spin-independent part 
as V'{r) = V{r) — Pave for the following reason: The 
quantity of V'{r) can be directly obtained from the BS 
amplitudes of IS'-charmonium states except the overall 
factor of I/mg. It ends up with less statistical uncer¬ 
tainties compared to the original potential V{r), whose 
estimation requires the subtraction of Pave- This is sim¬ 
ply because the value of Pave = -^ave ~ Smg receives 
somewhat large uncertainties, which arise mainly in the 
determination of mg through Eq. ( 8 ). Indeed, when we 
evaluate a difference between Vq and Pave directly from 
the fit with the Gornell functional form to the data of 
V'{r), this quantity shows much smaller error such as 
Vo — Pave = —0.146(13) GeV, in comparison to the val¬ 
ues of Pave = 0.508(69) GeV and mg = 1.789(34) GeV. 

We therefore adopt the energy-shifted potential of 
V'{r) to reduce uncertainties on the final result for energy 
eigenvalues, and then solve the following Schrddinger 
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equation for the reduced wave function usLjir) 


j usLj{r) = E'sLjUsLjir) 

I mq mqr^ j 

(15) 

where Vs^j(r) = VsLj{r) - E'ave and E'g^ j = Eslj - 
Eave with angular momentum quantum numbers {S, L 
and J). The interquark potentials Vgj^j{r), which may 
involve the spin-dependent interactions, clearly depend 
on the charmonium states labeled with specific S, L and 
J. The rest mass energy of the desired charmonium state 
is obtained simply by adding the energy eigenvalue of 
E'gLj to the spin-averaged IS'-charmonium mass of M^ve, 
which is evaluated by the standard lattice spectroscopy 
with high accuracy: Mslj = M^ve + E'gj^j = 2mQ + 
Eslj- 

The potential calculated from lattice QCD with the 
BS amplitude method are by definition discretized in 
space. In this context, instead of the continuum 
Schrodinger equation, we practically solve eigenvalue 
problems of finite-dimensional vector = u(nd) and 
finite-dimensional matrix [76] as 

^ ) Em.n^n — E Um- ( 16 ) 

n>0 


Note that a summation of n does not include n = 0 
since the reduced wave function is required to vanish at 
the origin. In the symmetric matrix H^ n for n,m > 0, 
diagonal and off-diagonal matrix elements are given by 


H-n. _ n — 


a^rriQ 

1 


L{L + 1) 


EnJzl^n — : 

a^iriQ 


+ V'{nd), (17) 
(18) 


and all other elements are zero. Here we omit the labels 
SLJ for clarity. 

In this work, we separately solve Eq. (16) in the direc¬ 
tions of vectors r which are multiples of (1, 0, 0), (1,1, 0) 
and (1,1,1). We prefer to use mainly on-axis data, which 
is expected to receive smallest discretization error and 
correction due to rotational symmetry breaking as stud¬ 
ied in Ref. [21], and take the largest difference between 
on-axis and off-axis results as the systematic error due 
to the choice of the r direction, while statistical errors 
are estimated by the jackknife method. Systematic un¬ 
certainties stemming from the choice of the fitting win¬ 
dow in the averaging process over the time-slice range are 
smaller than other errors. 

Alternatively, we may solve the continuum Schrodinger 
equation with the parameterized charmonium potential 


® Here, we assume that the reduced wave function vanishes at 
origin lim,— r(p'{r) = u{r) = 0. Indeed, if the potential satisfies 

r^V{r) 0, one can easily show the reduced wave fanction 

asymptotically behaves as u{r) 



FIG. 5: The energy-levels (dotted lines) and corresponding 
reduced wave functions u{r) (squares) of spin-averaged IS- 
charmonium states, obtained by solving the Schrodinger equa¬ 
tion with the lattice inputs. Only the central (average) val¬ 
ues are shown for both quantities. A horizontal line indi¬ 
cates the open-charm threshold. The spin-independent char¬ 
monium potential obtained from lattice QCD and its fitting 
result with the Cornell-plus-log functional form as a function 
of r are also overlaid in the same plot as circles and a shaded 
band, respectively. 


by empirical functional forms such as the Cornell form 
or the Cornell-plus-log form as discussed in Sec. IV. This 
procedure, however, yields large uncertainties in the low- 
lying energy levels, which highly depend on the choice of 
functional forms especially at short distances. To avoid 
such model dependence, we adopt the former strategy, 
which does not suffer from the sensitivity to the shape of 
the potential at short distances. 


B. Wave functions solving the Schrodinger 
equation 

In Fig 5, we plot energy-levels (dotted lines) and cor¬ 
responding reduced wave functions (curves with square 
symbols) up to the second excited state for the spin- 
averaged 5'-wave charmonium states, which are given by 
the charmonium potential with an expectation value of 
the spin operator being zero, {Sq ■ Sq) = 0. The spin- 
independent charmonium potential, which is calculated 
from lattice QCD, is also overlaid in the figure as circle 
symbols together with its fitting result using the Cornell- 
plus-log form (shaded band). 

We first carefully examine the energy eigenvalue E'^^^ 
of the spin-averaged IS” charmonium state whose mass 
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was used as input to calibrate the RHQ parameters for 
the charm quark. Recall that the value of i^ave is sup¬ 
posed to be zero because of its definition on the shifted 
energy E' introduced in Eq. (15). We consequently ob¬ 
tain i?ave = 0.2(1.3)(0.5) MeV, where the first error is 
statistical, the second error is systematic error due to ro¬ 
tational symmetry breaking. The obtained value is suf¬ 
ficient for satisfying the condition = 0 as a self- 

consistency check in our approach. We then conclude 
that the spin-averaged 15'-charmonium state can be well 
described by our charmonium potential given in the range 
r < 1.1 fm. 

The boundary condition implemented in the defini¬ 
tion of the Hamiltonian matrix defined in Eq. (16) en¬ 
forces the wave functions to vanish outside the interval 
r < dNs/2. Although our choice of Ns/2 x Ns/2 for 
the size of Hamiltonian matrix Hn^m is large enough for 
the IS'-charmonium states as discussed above, the higher- 
lying states that have more extended wave functions seem 
to suffer from the finite size effect caused by the bound¬ 
ary condition. Indeed, the resulting wave functions of 
the 2S and 3S charmonium states might be somewhat 
squeezed due to the smaller size of dNs/2. Therefore, 
these energy levels would be pushed down slightly due to 
the shrinkage of wave functions being less affected by the 
confining potential. As we mentioned above, the lattice 
data of the spin-independent potential becomes noisy in 
the range r > 1.1 fm, where signal-to-noise ratio of the 
BS wave function is poor, and also suffers from the fi¬ 
nite volume effect in lattice QCD simulations. In order 
to draw a firm conclusion for properties of higher-lying 
charmonium states without these effects, we clearly need 
to extend the calculation of the charmonium potential de¬ 
rived from the ground state wave function to the higher- 
lying states such as 2S and 3S states, which have more 
extended wave functions, using a sufficiently large lattice. 


C. Chromium mass spectrum 

We show the charmonium spectrum below 4200 MeV in 
Fig. 6. Theoretical spectra plotted as rectangular shaded 
boxes are given by solving the discrete nonrelativistic 
Schrodinger equations with the theoretical inputs. Verti¬ 
cal box length represents the level of uncertainty, which 
is given by adding statistical and systematic errors in 
quadrature. For the purpose of comparison, we plot both 
experimental values (horizontal lines) and results of the 
standard lattice spectroscopy (square symbols) together. 
The experimental values are taken from Particle Data 
Group [45]. At first glance, one can find that below the 
open charm threshold, our theoretical calculations from 
the NRp model with the lattice inputs excellently agree 
with not only the lattice spectroscopy, but also experi¬ 
ments. Especially a agreement between two lattice re¬ 
sults provides a strong check for the validity of our new 
method. All results including the lattice spectroscopy re¬ 
sults are also summarized together with the experimental 


TABLE VI: Masses of the charmonium states below 4200 
MeV are summarized in units of MeV. The labels of AVE 
and HYP in a column of “state” for 5'-wave states denotes the 
spin-averaged mass (Mi + SMa ) /4 and hyperfine splitting 
mass Msg^ — Mig^. Experimental data (denoted as Exp.) 
in the second column are taken from Particle Data Group, 
rounded to 1 MeV [45]. There are two kinds of lattice QCD 
results tabulated in the third and fourth columns. One is ob¬ 
tained by the standard lattice spectroscopy, while another is 
evaluated by solving the Schrodinger equation with the char¬ 
monium potential determined from lattice QCD. For the lat¬ 
ter, the first error is statistical and the second error systematic 
as described in text. The spin-weighted average mass (de¬ 
noted as ^[L]j) are also included for spin triplet states ^[E\j. 
The last column shows the results from a NRp model [4]. 


state Exp. 

Lattice QCD 

NRp model [4] 



spectroscopy BS amplitude 


Vc (FSo) 

2981 

2985(1) 

2985(2)(1) 

2982 

j/v> (i^Si) 

3097 

3099(1) 

3099(2)(1) 

3090 

AVE 

3068 

3070(9) 

3070(2)(1) 

3063 

HYP 

116 

114(1) 

113(1)(0) 

108 

Vc (2^So) 

3639 


3612(9)(7) 

3630 

(2®Si) 

3686 


3653(12)(5) 

3672 

AVE 

3674 


3643(11)(5) 

3662 

HYP 

47 


41(6)(3) 

42 

Vo [3^So) 



4074(20)(70) 

4043 

V' (3®Si) 

4039 


4099(24)(98) 

4072 

AVE 



4092(22)(91) 

4065 

HYP 



25(15)(28) 

29 

ho (HPi) 

3525 

3506(6) 

3496(7)(19) 

3516 

(DPj) 

3525 


3503(7)(10) 

3524 

Xco (l^Fo) 

3415 

3393(6) 


3424 

Xci (I'Ri) 

3511 

3485(6) 


3505 

Xc2 (1"P2) 

3556 



3556 

ho (2iPi) 



3927(16)(34) 

3934 

(23Pj) 



3916(19)(31) 

3943 

XcO (2®Po) 

3918 



3852 

Xci (2®Pi) 




3925 

Xc2 { 2 ^P 2 ) 

3927 



3972 

Vc 2 (1^1)2) 



3783(12)(4) 

3799 

(DP,/) 



3774(13)(2) 

3800 

(i®Pi) 

3773 



3785 

V' (1®P2) 




3800 

(I'Ps) 




3806 

Vo 2 (2^02) 



4221(21)(72) 

4158 

(23p,/) 



4193(25)(88) 

4159 

(2®Pi) 

4153 



4142 

(2^02) 




4158 

(2"P3) 




4167 


values in Table VI. 

In this study, we have succeeded in extracting only 
the spin-spin potential among the spin-dependent parts 
of the interquark potential. Thus at this stage we cannot 
predict the spin-orbit splitting which is led by the tensor 
and spin-orbit terms of the spin-dependent potential. In 
other words, we can compute only the spin-averaged mass 
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FIG. 6: Mass spectrum of charmonium states below and near the open-charm threshold. The vertical scale is in units of MeV. 
Labels of are displayed in the lower (upper) horizontal axis. Rectangular shaded boxes indicate predictions from 

the NRp model with purely theoretical inputs based on lattice QCD and their errors which are the sum of the statistical and 
systematic added in quadrature. Solid lines indicate experimental values of well established charmonium states, while square 
symbols represent results of the standard lattice spectroscopy. A horizontal solid line shows the open-charm threshold. A symbol 
of ^Pj denotes the spin-weighted average of spin-triplet ^Pj states whose mass is given by -|- )/9- 


for higher-wave charmonium states like P-wave charmo¬ 
nium Xcj state. 

The mass splitting between the radial excitations and 
the ground state also provides an important validity 
check on our new approach. Fig. 7 shows several mass 
splittings theoretically predicted by the hybrid approach 
in comparison to physical values of the corresponding 
splittings. In the top left panel of Fig. 7, the radial 
excitation mass splitting of the spin-averaged 15' and 
25 states is Mave(25) - Mave(15) = 573(10)(5) MeV, 
of which value is slightly smaller than the experimental 
value of 606(1) MeV [45]. This deviation (~ 30 MeV) 
from the experiment can be attributed to the finite vol¬ 
ume effect, which is caused by the fact that more ex¬ 
tended 25 state than 15 state is forced to fit in the spa¬ 
tial volume ~ (3 fm)^. Note that the S-D mixing due 
to the tensor force is simply ignored in our calculations. 
Thus, the spin-spin potential solely gives the mass split¬ 
ting in hyperfine multiplets. The top right panel of Fig. 7 
shows the hyperfine-splitting energy for the 15 charmo¬ 
nium states: Mis^yp = ~ ^ri^- H is found that 

there is a good agreement between two lattice results; 


113.4(9)(1) GeV from the NRp model with lattice inputs 
and 113.8(8) GeV from the standard lattice spectroscopy. 
This simply suggests that there is no additional uncer¬ 
tainties induced by both determining the charmonium 
potential and the charm quark mass by the BS ampli¬ 
tude method and solving the nonrelativistic Schrodinger 
equation with them. 

We remark that the values of the hyperfine-splitting 
are slightly smaller than the experimental one = 

116.6(1.2) GeV. This would be simply due to insufficient 
calibration of the RHQ parameters and also other pos¬ 
sible systematic uncertainties including the remnant dis¬ 
cretization artifact. On the other hand, the hyperfine 
splitting energy for the 25 charmonium states, which is 
plotted in the bottom left panel of Fig. 7, shows that 
the value of M^/(2S) — M^X 2 S) = 41(6)(3) MeV obtained 
from the hybrid approach is roughly consistent with the 
experimental value 47(1) MeV, within its error range. 

The bottom right panel of Fig. 7 shows the IP hy¬ 
perfine mass splitting which is given by an energy dif¬ 
ference between between the and spin-averaged Xcj 
states: Mip^hyp = ~ Experimentally, the 
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FIG. 7: Mass splittings of states lying below the open 
charm threshold in units of MeV, compared to the phys¬ 
ical mass splitting. Upper panels show the mass splitting 
Mave(2S') — Mave(lS) between the spin-averaged IS- and 2S- 
states (left), and the hyperfine mass splittings Mj/^ — Mn^ 
between IS-states (right), while lower panels show mass split¬ 
ting of M^( 2 S) ~ Mn^( 2 S) between 2S-states (left) and mass 
splitting of — Mh^ between IP-states (right). In each 

plot, cross and circle symbols indicate the experimental data 
and theoretical results obtained from the NRp model with the 
lattice inputs, respectively. The quoted errors indicate the 
sum of the statistical and systematic errors added in quadra¬ 
ture. Dashed lines represent the central value of the exper¬ 
iment. Only for the hyperfine mass splitting between 15- 
states, the results of the lattice spectroscopy is shown as a 
square. 


value of is known to be zero with high accu¬ 
racy as p = 0.02(23) MeV [62, 63]. The hybrid 

approach yields a small splitting energy as Mip^hyp = 
7.2(1.6)(9.3) MeV, which is consistent with the zero value 
within a large error. Of course, however, the spin-spin 
charmonium potential determined in the BS amplitude 
method is not still enough to describe the tiny IP hyper¬ 


fine splitting measured in experiment. As we mentioned 
in the previous subsection, a finite-range spin-spin po¬ 
tential gives a nonzero value of hyperfine mass splitting 
even in the case of higher-wave states such as P-wave 
state, while zero hyperfine splitting measured in exper¬ 
iments is easily reproduced by the point-like spin-spin 
potential widely adopted in phenomenological quark po¬ 
tential models. Here, we stress that the spin-spin poten¬ 
tial from the BS amplitude method is finite-range and 
therefore the value of is highly sensitive to both 

shapes of the spin-spin potential and wave functions of 
P-wave states. According to our systematic study of 
the BS amplitude method performed in quenched lattice 
QCD [21], the spin-spin potential receives large uncer¬ 
tainties due to the discretization artifacts more than the 
spin-independent central potential. To make a firm con¬ 
clusion, it is necessary to perform the present calculation 
on the finer lattice. 

Our theoretical calculations for the charmonium mass 
spectrum below the open-charm threshold are basically 
in good agreement with the experimental measurements. 
The point we wish to emphasize here is that our novel ap¬ 
proach has no free parameters in solving the Schrodinger 
equation in contrast to the phenomenological NRp mod¬ 
els. All of the parameters are fixed by lattice QCD simu¬ 
lations, where three light hadron masses {e.g. pion, kaon 
and fl baryon) are used for setting the lattice spacing a 
and hopping parameters of the light and strange quarks 
{i.e. the light and strange quark masses). In this study, 
the charm quark was treated in the quenched approxima¬ 
tion. Then the experimental values of rjc and J/^/> char¬ 
monium masses are used to determine the charm quark 
parameters appeared in the RHQ action. In this sense, 
the hybrid approach proposed here is distinctly different 
from existing calculations in the phenomenological quark 
potential models. 

Let us now attempt to straightforwardly extend the hy¬ 
brid approach to above the open-charm threshold. Only 
the spin-averaged mass is co nsidered for the P and D 
spin-triplet states: M{n^Pj) = (M„3 Po + 3M„3pj -I- 
5Mn3pJ/9 and MfnfDj) = (3M„3p)^ -|- + 

7M„ In order to provide mass splittings among 

these spin-triplet states, the tensor and spin-orbit poten¬ 
tials are inevitably required. Since, in this paper, we 
succeeded in extracting the spin-spin potential solely for 
the spin-dependent potentials, we should focus on the 
spin-averaged masses. 

First of all, one can observe that the values obtained 
from the hybrid approach above the open-charm thresh¬ 
old fairly agree with the existing experimental data, al¬ 
though errors are relatively large as shown in Table VI. 
We, however, are not in a position to give a realistic de¬ 
scription to the higher-lying charmonium states, which 
are located above the open-charm threshold. This is be¬ 
cause there are the following remarks in our calculations 
including the higher-lying charmonium states. 

1. The higher-lying charmonium states significantly 
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suffer from systematic uncertainties, which are 
mainly due to the less knowledge of the long-range 
part of the spin-independent potential. We have no 
reasonable data for the charmonium potential at 
longer distances than about 1.1 fm since the wave 
function of the IS' ground-state possesses the highly 
localized nature. Therefore we need to calculate the 
potential form the higher-lying charmonium states. 
Alternatively, we simply extrapolate the long-range 
behavior of the potential outside the region, where 
the charmonium potential is really determined from 
the localized wave function. In the latter case, the 
higher-lying spectrum of the charmonium is more 
sensitive to the choice of the adopted functional 
form in the fitting procedure. 

2. The possible mass shift due to mixing the QQ states 
with DD continuum is completely neglected in this 
study. One may expect that the NRp models with¬ 
out such mixing works well to describe the low- 
lying charmonium systems far below its threshold. 
On the other hand, such coupled channel effects 
might not be negligible near and above the thresh¬ 
old and then the potential description may lose 
the accuracy of theoretical prediction, though the 
naive treatment of the NRp models even for higher- 
lying charmonium states was phenomenologically 
successful despite the absence of coupled channel 
effects [3, 4] 

3. For the higher-lying excitations of the spin-1 char¬ 

monium state, the S-D mixing becomes severe 
since the level spacings between (n -|- and 

n^Di get narrower [77]. However, S-D mixing ef¬ 
fects on J/V’, ■*/'(25') and states are not taken 

into account in the present calculation since the 
tensor term in the spin-dependent potentials is not 
determined in this study. Similarly, the mass es¬ 
timations of Xcji^P) and iplnD) tabulated in Ta¬ 
ble VI are calculated without consideration of pos¬ 
sible partial-wave mixings such as S-D^ F-P and 
D-G mixings. 

To calculate the BS wave function of IP-states, better 
source operators with respect to odd-parity wave func¬ 
tion [78] are required. Meanwhile some extension of the 
variational method [79, 80] to the four-point correlation 
functions is necessary for extracting the BS wave function 
of the radial excitation of the 5'-wave states. These new 
calculations can give more realistic prediction especially 
to the higher-lying charmonia. The former provides in¬ 
formation of the spin-orbit and also tensor potentials [78]. 
The latter can provide not only the tensor potential, but 
also the mixing angle between 2^S'i and l^Di states in 
the same way as the nuclear force [25] . Furthermore more 
data points of the charmonium potential at large dis¬ 
tances can be accessible from such excited states of the 
charmonium, which have more extended wave function 
than that of IS” ground states. Such kind of studies is 
now under way [81]. 


TABLE VII: Masses of low-lying Ds meson states, the spin- 
averaged mass and hyperfine splitting energy of IS" charmo¬ 
nium states. The colums have the same meaning as in Ta¬ 
ble VI. Results are given in units of GeV. 


state 

{jn 

r 

fit range 

mass [GeV] 

xVd.o.f. 

Ds 

(0-) 

75 

[30:47] 

1.9780(12) 

1.08 

Dl 

(1-) 

H 

[30:47] 

2.1230(42) 

0.61 

Mave(lS) 

— 


— 

2.0865(33) 

— 

Ahyp(15) 

— 


— 

0.1461(37) 

— 

Dio 

16^ 

1 

[14:26] 

2.3536(77) 

1.45 

Dsi 

(1+) 

757i 

[14:26] 

2.4689(83) 

1.14 

Dsi 

(1+) 

7i7j 

[14:22] 

2.4893(87) 

1.20 


VI. APPLICATION TO HEAVY-LIGHT 
SYSTEM 

In the charmonium (heavy-heavy) system, the spec¬ 
trum below the open charm threshold are well described 
by potential description with our charmonium poten¬ 
tial including the spin-spin interaction, which was de¬ 
termined by the BS amplitude method in dynamical lat¬ 
tice QCD simulations. In this section, we apply the new 
method to the Dg heavy-light meson system, which rep¬ 
resents the case of mesons with non-degenerate quark 
masses. Apart from the phenomenological interest, we 
also would like to examine the validity range of the new 
method in terms of the size of quark kinetic mass. 

A. Lattice setup for charmed-strange mesons 

The numerical setup for the charmed-strange {Ds) 
mesons system is basically same as for the calculation 
of the charmonium system. For the charm quark, we 
use the RHQ action with the parameters calibrated by 
IS'-charmonium states. In addition to the computation 
of the charm quark propagator, another fermion matrix 
inversion for a strange quark is required to compute the 
Ds-meson correlation functions. The non-perturbatively 
O(a)-improved Wilson quark action {csw = 1-715) is 
used for the strange quark. A hopping parameter of the 
strange quark is chosen to be Ks = 0.13640, which is the 
same as the sea strange quark used in gauge held gen¬ 
eration. Simultaneously, ss-mesons are supplementarily 
calculated, and we use sAmeson data for a consistency 
check on the kinetic masses of the strange quark deter¬ 
mined through ss and cs systems as we will discuss later. 

Fig. 8 shows the effective mass plots for S and P-wave 
Ds-meson states. The Pg-meson masses are determined 
by constant hts to the plateaus observed in the effec¬ 
tive mass plots with covariance matrices accounting for 
the data correlation among different time slices. Re¬ 
sults of the Pg-meson masses together with ht ranges 
used in the hts and the values of x^/d.o.f are summa¬ 
rized in Table VII. The quoted errors represent only the 
statistical errors given by the jackknife analysis. The 
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FIG. 8: Effective mass plots for low-lying S and P-wave Ds 
meson states. P-wave states {Ds and D*) and P-wsve states 
(Z)so(2317), Psi(2460) and Psi(2536)) are shown in the up¬ 
per and lower panels, respectively. Each Da meson state is 
specified in legend. Horizontal lines and shaded bands denote 
fit results with statistical errors estimated by the jackknife 
method and their fit range. 


RHQ action for the charm quark works well even for the 
low-lying Zlg-mesons. The spin-averaged and hyperfine 
splitting Hg-meson masses, = 2.0865(33) GeV and 
Mj(yp = 0.1461(37) GeV, are obtained from the standard 
lattice spectroscopy. Although the simulated strange 
quarks are slightly off the physical point, these results 
are quite close to the experimental data of Mf(l|’(lS') = 
2.07635(27) GeV and M“P(1S') = 0.1438(4) GeV. The 
deviations from the experimental results are within about 
0.5%. Furthermore, results of P-wave Pg-meson states 
from the lattice spectroscopy marginally reproduce the 
experimental data. Similar results are reported by the 
PACS-GS collaboration using 2-1-1 flavor dynamical 
gauge configurations generated with the physical strange 
quark [44, 82]. 

The two-point correlation functions of both pseu¬ 
doscalar and vector ss-mesons, i.e. r7ss(0“) and (()(1“) 
mesons, are also calculated in this study. We ob¬ 
tain results of = 0.7699(9) GeV and = 

1.0827(68) GeV. The similar values are reported in 
Ref. [22]. The fit range was chosen to be 24 < t < 
39 for both states. The the (j) meson mass is some¬ 
what heavier than the experimental values of = 



FIG. 9: The reduced BS wave function u{r) = r(j){r) for the 
IS' vector cc (circles), cs (squares) and ss (triangles) states, 
as a function of spatial distance r. They are normalized as 

E<^(®) = 1- 

1.019455(20) GeV. It should be attributed to the fact 
that the simulated strange quarks are slightly off the 
physical point. Although the systematic uncertainty due 
to slightly heavier strange quark mass is expected to be 
extremely small in the charmonium spectrum, we should 
take into account some corrections for the Pg-meson 
spectrum [46]. 


B. BS wave function 

In Fig. 9, we show the reduced BS wave functions for 
the IS” vector cc, cs and ss states corresponding to J/i), 
D* and (j) mesons, respectively. It is found that the D* 
wave function is spatially extended to at least the half 
of the spatial extent of lattice volume (Afga/2 ~ 1.5 fm). 
Although the amplitude of the wave function of the P* 
meson is considerably small at r ~ 1.5 fm, it still seems 
to remain non-zero values in the range of r > 1.5 fm, 
where only off-axis data points are available. The wrap 
round effect would cause the rotational symmetry break¬ 
ing at longer distances. Therefore, in the Pg system, 
the interquark potential could be more affected by the 
finite volume effect than the charmonium system. In the 
case of the ss system, which is more spatially extended 
than the cs system as shown in Fig. 9, this problem could 
become more severe. 


C. quark kinetic mass 

Fig. 10 illustrates the determination of quark kinetic 
mass from the cs and ss meson systems in the BS ampli¬ 
tude method. A quantity iriQq is defined as twice the re¬ 
duced mass of the Pg(cs) system: toq, = 2mQmq/{mQ + 
niq), while corresponds to the strange quark mass. 

As shown in the upper panel of Fig. 10, we fit the 
data points of the cs meson system at relatively large 
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FIG. 10: The determination of the reduced mass mqq 

from the Ds{cs) system (upper) and strange quark mass niq 
from the ss system (lower) in the BS amplitude method. 
We obtain the quark kinetic masses of rriQq and niq from 
the asymptotic behavior of the right-hand side of Eq. (8) in 
long-distance region. Solid lines with shaded bands repre¬ 
sent the fitting results and ht ranges with the statistical er¬ 
ror estimated by the jackknife method. In the lower plot, 
a horizontal shaded band indicates a kinetic mass of the 
strange quark, which is independently evaluated by the re¬ 
lation TUq = mQqmQ/(2mQ — ruQq) with the values of mqq 
and mq from the cs and cc systems. 


distances, where the reasonable plateau is found in the 
region of r > 0.7 fm. We then obtain the value of 
i^Qq = 0.959(45) (34) (36) GeV. The first error is sta¬ 
tistical, and the second and third ones are systematic 
uncertainties due to the choice of data points taken 
from three directions and a variation of tmin, respec¬ 
tively. The strange quark mass nig can be evaluated 
by two data sets of mg^ and mg through the relation 
niq = niQqmQ/{2niQ — niQg). The value of mg corre¬ 
sponds to the charm mass, which was already determined 
in the previous section. When combined with results ob¬ 
tained from the Ds{cs) and charmonium (cc) system, we 
obtain the value of niq = 656(41) MeV for the strange 
quark mass. Quoted error is statistical one, which was 
determined by the jackknife method. 

Independently, the strange quark mass niq can be de¬ 
termined through the ss system as depicted by circle 
symbols in the lower panel of Fig. 10. Similar to the 


TABLE VIII: Charm and strange quark masses, which are 
determined from the Coulomb-gauge quark-antiquark BS am¬ 
plitude and the Landau-gauge quark propagator, are summa¬ 
rized in units of GeV. 



BS amplitude 

quark propagator 

flavor 

QQ or qq Qq 

Landau gauge 

charm 

1.784(23) — 

1.776(8) 

strange 

0.554(19) 0.656(41) 

0.643(5) 


upper hgure, the reasonable plateau is found in the re¬ 
gion of r > 0.7 fm. We compute a weighted average of 
the data points in the range of 8 < r/a < 7^/3 with 
a covariance matrix accounting for the correlation, and 
then obtain the value of niq = 554(19) (6) (8) MeV, which 
is close to a typical value of constitute strange quark 
mass (~ Mrj,/2 « 500 MeV) adopted in SU{6) quark 
models [1]. The meaning of the three quoted errors is 
explained above. 

For comparison, the previously estimated value of 
niq = 656(41) MeV from the data sets of mg, and mg is 
also displayed by a horizontal shaded band in the lower 
figure. There is 2cr discrepancy between this band and 
the plateau behavior of —(V^(/)v/^^>v —V^()>ps/(/>ps)/£'hyp 
for the ss system at large distances. Although this dis¬ 
crepancy may imply that nonrelativistic treatment is no 
longer valid for the heavy-light system, we would like to 
remind of the fact that the BS wave function of the ss- 
meson system at large distances is likely affected by the 
finite volume effect as discussed previously. 

The following discussion shows that the above specu¬ 
lation is likely to be true. The strange quark mass deter¬ 
mined from the BS wave function of the ss-meson states 
is indeed underestimated compared to a “pole mass” de¬ 
termined from the effective mass of gauge-variant quark 
two-point correlator in the Landau gauge, while two dif¬ 
ferent calculations show remarkable consistency for the 
charm quark as summarized in Table VIII. Fig. 11 shows 
effective mass plots and comparisons with ones obtained 
within the BS amplitude method. The Landau-gauge 
pole mass is an alternative way of measuring the quark 
mass in lattice QCD. For details of how to calculate it, 
see Ref. [83]. 

If one may choose the cs result rather than the ss result 
for the BS amplitude method, the strange quark masses 
from two estimation methods become consistent again. 
Although the physics behind the consistency discussed 
here is beyond the scope of this paper, we may simply 
conclude that the discrepancy between results from the 
cs and ss mesons is mainly attributed to the finite volume 
effect on the ss wave function. 


D. charmed-strange potential 

Fig 12 shows results of the spin-independent and spin- 
spin interquark potentials obtained from the Ds and D* 
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time slice t/a 



time slice t/a 


FIG. 11: Effective masses of gauge-variant quark two-point 
correlator in the Landau gauge for charm (upper panel) and 
strange (lower panel). Solid lines indicate fit results for “pole” 
masses of charm and strange quarks and shaded bands display 
the fitting ranges and one standard deviations estimated by 
the jackknife method. In each panel, a wider and horizontal 
shaded band indicates a kinetic mass evaluated from the cs 
and/or cc systems within the BS amplitude method. 

meson states (hereafter called cs potential) in the dy¬ 
namical QCD simulation. For purpose of comparison, 
the charmonium potentials are also displayed. 

At hrst glance, a shape of the cs potential is basically 
similar to that of the charmonium potential, so that we 
similarly adopt the Cornell potential form for the spin- 
independent cs potential and also exponential (Yukawa) 
form for the spin-spin cs potential as is the case in the 
charmonium potential. We obtain the Cornell parame¬ 
ters of the cs potential as A = I.30(8)(22)(21)(21) and 
i/(T = 324(16) (34) (26) (4) MeV with a reasonably small 
value of x^/d.o.f. « 1.9. The first error is statistical and 
the second, third and forth ones are systematic uncer¬ 
tainties due to the choice of data points taken from three 
directions, and variations of tmin and Tmin, respectively. 

The appropriate fitting range was determined to min¬ 
imize a x^/d.o.f. value taking into account the data cor¬ 
relation among different spatial distances r. We choose 
the ht range of [rmin/o : ?’max/a] = [4 : 7-\/3], which 
corresponds to the same range in the case of the char¬ 
monium potential. Although the string tension has weak 
dependence on quark kinetic mass [21], the Coulomb co¬ 
efficient of the cs potential significantly grows in com- 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 


r [fm] 

FIG. 12: The spin-independent and spin-spin interquark po¬ 
tentials for the cc (circles) and cs (squares) systems, calcu¬ 
lated from the BS wave functions in dynamical lattice QCD 
simulations with almost physical quark masses. In the upper 
panel, the spin-independent parts of both the charmonium 
and cs potentials are plotted. For clarity of the figure, the 
constant energy shift J5ave, which is given by the spin-averaged 
mass of 15 states, is not subtracted. Solid and dashed curves 
represent the fit results with the Cornell parametrization. The 
shaded bands show statistical uncertainties in the fitting pro¬ 
cedure where the jackknife method is employed. In the lower 
panel, we show the spin-spin potential Vs{r). The exponen¬ 
tial form is used for fitting the resultant spin-spin potentials 
for the cc and cs systems. The inset shows a magnified view. 


parison with that of the charmonium potential. For the 
spin-spin potential, we obtain a = 3.79(36) GeV and 
P = 2.89(9) GeV (a = 1.48(14) and P = 1.97(9) GeV) 
from the exponential (Yukawa) form fit with y^/d.o.f. « 
1.49 (x^/d.o.f. « 1.86). We quote only the statistical 
errors, which are determined by the jackknife method. 
We find that the size of the finite-range of the spin-spin 
potential for the cs system is almost consistent with the 
one obtained from the charmonium spin-spin potential 
within statistical uncertainties. 
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FIG. 13: Mass spectrum of the charmed-strange mesons around the DK threshold. The vertical scale is in units of MeV. Labels 
of j (J^) are displayed in lower (upper) horizontal axis. Solid lines indicate experimental values of well established Dg 

meson states, while square symbols represent results of the standard lattice spectroscopy. Rectangnlar shaded boxes indicate 
predictions from the NRp models with purely theoretical inputs based on lattice QCD and their errors which are the sum of 
the statistical and systematic errors added in quadrature. A horizontal solid line shows the DK threshold. A symbol of ^Pj 
denotes the spin-weighted average of spin-triplet ^Pj states. 


E. charmed-strange meson mass spectrum 

Using the cs potential calculated from lattice QCD, we 
obtain the spectrum of the charmed-strange mesons in 
the same footing as the charmonium spectrum discussed 
in Sec. V C. The resulting Dg meson spectrum together 
with the experimental values and the results from the 
standard lattice spectroscopy are summarized in Table IX 
and also Fig. 13. 

The results of IS' states in the Dg-meson family by 
solving the discrete nonrelativistic Schrodinger equation 
with lattice inputs show a good agreement with both 
experiments and standard lattice calculations below the 
DK threshold. Above the DK threshold, the cs poten¬ 
tial obtained from the BS amplitude method can properly 
reproduce the mass ordering of the Dg mesons, while ab¬ 
solute values of the masses are systematically larger than 
the experimental values; for example, the corresponding 
Dgi(2536) state is overestimated by about 90 MeV, and 
then observe a large systematic discrepancy between two 
kinds of lattice QCD results. 

Recall that the results obtained from the standard lat¬ 
tice spectroscopy near the DK threshold are also slightly 


deviated from the experimental values. Although this im¬ 
plies that the observed discrepancies among them are not 
solely attributed to the BS amplitude approach, there are 
three possible sources of above mentioned discrepancies 
in our method. The first is, as we have noted repeatedly, 
associated with uncertainties in the long-range part of the 
spin-independent potential. The fact that the Dg state 
has the wider wave function leads to the large finite vol¬ 
ume effect on the interquark potential determined within 
the BS amplitude method at long distances. Further¬ 
more, information of the cs potential is accessible only 
within the “size” of the S'-wave Dg mesons. We therefore 
force the wave function to be zero at the spatial bound¬ 
ary in the process of solving the discrete nonrelativistic 
Schrodinger equation with present limited data of the cs 
potential. 

The second possibility is that there are the DK and 
D*K threshold effects. Since Dgo(2310) and Dgi(2460) 
are located near the DK and D*K thresholds respec¬ 
tively. Therefore, the coupling of these states to the 
DK and D*K two-hadron states could not be negligi¬ 
ble [85, 86]. Such channel couplings may cause level re¬ 
pulsion and thus mass shift by the threshold effect might 
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TABLE IX: Masses of low-lying Ds mesons are summarized 
in units of MeV. The labels of AVE and HYP in a column 
of “state” for S-wave states denotes the spin-averaged mass 
(Ml + 3 M 3 Sj )/4 and hyperfine splitting mass Ms— Mi . 
Experimental values in the second column are taken from 
Particle Data Group, rounded to 1 MeV [45]. There are 
two kinds of lattice QCD results tabulated in the third 
and fourth columns. One is obtained by the standard lat¬ 
tice spectroscopy, while another is evaluated by solving the 
Schrodinger equation with the cs potential determined from 
lattice QCD. For the latter, the first error is statistical and 
the second error systematic as described in text. The spin- 
weighted average mass (denoted as ®[L]j) are also included 
for spin triplet states ^[L]j. The last column gives the results 
from a NRp model [84] . 


state 

Exp. 

Lattice QCD (This work) 

NRp model [84] 



spectroscopy BS amplitude 


Ds (Us'd) 

1968 

1978(1) 

2000(10)(3) 

1963 

D: (l^^i) 

2112 

2123(4) 

2138(8)(3) 

2099 

AVE 

2076 

2087(3) 

2103(8)(3) 

2065 

HYP 

144 

146(4) 

138(5)(1) 

136 

Ds {2^So) 



2766(38)(50) 


D*s {2^Si) 

2709 


2857(42) (80) 


AVE 



2834(40)(73) 


HYP 



92(9)(30) 


Dsi (UPi) 

2535 

2489(9) 

2623(30)(32) 

2527 

Dsj (DPj) 

2506 


2629(29)(32) 

2532 

D:o (I'Bo) 

2318 

2354(8) 


2446 

Dsi (l^Pi) 

2460 

2469(8) 


2515 

D:^ (UPz) 2572 



2561 


be happened [87-89]. 

Finally, it is worth pointing out that we assumed that 
there is no S-D mixing due to the tensor force because 
of large energy gap between and ID states. Strictly 
speaking, however, the D* meson which is specified by 
the quantum numbers = 1“ is not purely composed 
of the ^Si wave function. This approximation could in¬ 
troduce a small correction to the intermediate and short- 
range parts of the cs potential calculated in the BS am¬ 
plitude method. 

Therefore, we need further development of our ap¬ 
proach to take into account both the coupled channel ef¬ 
fect and S-D mixing. A simulation with sufficiently large 
volume is also required for precise prediction of masses 
of the Z?s-meson states near and above the DK {D*K) 
threshold within our approach. 


VII. SUMMARY 

We have calculated the interquark potentials for both 
the charmonium (cc) and charmed-strange (cs) mesons 
at almost physical point. The interquark potential with 
finite quark masses are defined through the equal-time 
Bethe-Salpeter wave functions of the pseudoscalar and 
vector mesons. Our simulations have been done in the 


vicinity of the physical light quark masses, which corre¬ 
sponds to Mjr « 156 MeV, using the PACS-CS Iwasaki 
gauge configurations with 2-I-I flavors of dynamical clover 
light quarks. We use the relativistic charm quark tuned 
to reproduce the experimental values of the ijc and J/tp 
masses. 

We first investigated the charmonium potential. The 
resulting spin-independent potential has the Coulomb- 
plus-linear form, and their parameters are close to the 
values used in the phenomenological NRp models. The 
string breaking due to the presence of dynamical sea 
quarks is not apparently observed. The spin-spin po¬ 
tential obtained from the dynamical simulations exhibits 
the finite-range repulsive interaction. Its shape is quite 
different from a repulsive ^-function potential induced 
by the one-gluon exchange, which is often adopted in the 
NRp models. 

Our ultimate goal is to reveal the mystery behind 
rich structures recently observed in the heavy-heavy 
and heavy-light systems including the newly discovered 
charmonium-like mesons. As a first step, we calculated 
the charmonium mass spectrum by solving nonrelativis- 
tic Schrodinger equation with purely theoretical inputs of 
the spin-independent and spin-spin potentials, and also 
the quark kinetic mass. To avoid any model dependence 
from fitting, we practically solve the discrete Schrodinger 
equation in finite volume with Dirichlet boundary condi¬ 
tion, and thus can handle direct lattice data of the char¬ 
monium potential without any parameterization. 

We found an excellent agreement of low-lying charmo¬ 
nium masses between our results and the experimental 
data. We here emphasize that our novel approach has 
no free parameters in solving the Schrodinger equation 
in contrast to the phenomenological NRp models. In our 
calculations, three light hadron masses {e.g. the pion, 
kaon and Q baryon are chosen in the PACS-CS collab¬ 
oration) and two charmonium masses (the rjc and J/V') 
are used for fixing the lattice spacing and quark mass 
parameters in the lattice QCD action including the RHQ 
parameters for the charm quarks. 

In order to precisely predict the mass spectrum above 
the open charm threshold, we should take into account 
both coupled-channel effect with the DD continuum and 
S-D mixing due to the presence of the tensor force. 
In this study, we simply ignore these effects and then 
apply the hybrid approach to higher-lying charmonium 
states above the open-charm threshold. We found that 
the theoretical predictions of the NRp model calculation 
with lattice inputs are remarkably consistent with well- 
established experimental data for the conventional char¬ 
monium states. 

For an application, we straightforwardly extend our 
method to calculate the charmed-strange meson system, 
which represents the case of mesons with non-degenerate 
quark masses and also heavy-light system. A shape of the 
interquark cs potential in the Ds-meson system is basi¬ 
cally similar to that of the charmonium potential. Using 
the resulting cs potential as theoretical inputs, we obtain 
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the spectrum of the charmed-strange mesons. Below the 
DK threshold, our new method works well in spite of 
the fact that the mesons contain a strange quark or 
strange anti-quark. Although above the DK threshold 
our cs potential can reproduce the mass ordering of the 
Ds mesons, absolute values of the masses are consistent 
with neither the experimental values nor the results from 
the standard lattice spectroscopy. 

Although it is difficult to draw any firm conclusion 
at this stage, the discrepancies observed above the DK 
threshold suggest that all of coupled-channel effect, S-D 
mixing and higher-order relativistic corrections, that are 
omitted in this study, are possibly important to under¬ 
stand the properties of the heavy-light mesons near the 
threshold. To disentangle these effects, simulations in 
a larger lattice is necessary. Conversely, one might say 
that by taking into account these effects properly, our ap¬ 
proach can give a systematic way to examine the validity 
of the potential description even for the charmed-strange 
system. Such full analysis finally sheds light on the de¬ 
tailed properties of the D^ mesons. 

At least, in this study, the charmonium and charmed- 
strange potentials obtained from the BS amplitude 
method have been succeed in reproducing the low-lying 
masses below the open charm threshold and DK thresh¬ 
old respectively. Furthermore we showed that our new 
analysis can potentially shed light on the detailed prop¬ 
erties of the heavy quarkonium system. While only en¬ 
ergy eigenvalues are evaluated from temporal informa¬ 
tion of meson correlation functions in the standard lat¬ 
tice spectroscopy, the new method takes an advantage of 
full spatial information together with temporal informa¬ 
tion. The BS wave functions can be identified with the 
eigenstates of the Hamiltonian. Hence, without knowing 
the details of an explicit form of the Hamiltonian, lots of 
physical quantities could be calculated directly by the BS 
wave functions as studied in the NRp models. For exam¬ 
ples, El and Ml radiative partial widths are supposed 


to be evaluated with the BS wave functions of the char¬ 
monium states. Such information is important to reveal 
the structure of hadrons. 

To derive a complete nonrelativistic Hamiltonian of the 
heavy-heavy and heavy-light systems from lattice QCD, 
we must calculate all spin-dependent terms (spin-spin, 
tensor and spin-orbit forces), which are required for more 
realistic predictions for the higher-lying states. We now 
develop the BS amplitude method to calculate the BS 
wave functions of P-wave mesons, which provides infor¬ 
mation of the spin-orbit and also tensor potentials. 

Once all spin-dependent terms of the interquark poten¬ 
tial are determined and also all systematic uncertainties 
are well understood, we will gain new and valuable in¬ 
sight on the mesons newly discovered in the heavy-heavy 
and heavy-light systems. It is also important to examine 
the validity of the potential description with the BS am¬ 
plitude method from the viewpoint of the ^-expansion 
for a nonlocal and energy-independent interquark po¬ 
tential U originally appeared in Eq. (4). For this pur¬ 
pose, we would like to examine whether the same in¬ 
terquark potential is obtained from the BS wave function 
of the radial excitation of the S'-wave states. All of the 
above-mentioned extensions of the new method are now 
in progress [81]. 
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